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Abstract 

Given  an  axially  symmetric  blunt  body  in  a  steady  super- 
sonic flow  of  compressible  fluid,  the  problem  is  to  study  the 
flow  between  the  body  and  the  detached  shock  wave  formed  in 
front  of  it.   Since  the  location  of  the  shock  is  unknown,  the 
inverse  problem  is  solved,  i.e.,  given  an  analytic  shock  shape, 
we  compute  the  flow  and  generating  body. 

Hyperbolicity  of  the  partial  differential  equations  of 
motion  is  achieved  in  the  subsonic  region  by  means  of  a  complex 
extension  of  the  initial  data.   The  resulting  equations  are  then 
solved  both  by  a  three-dimensional  method  and  by  the  method  of 
characteristics.   The  method  of  characteristics  is  found  to 
yield  more  information  about  the  flow. 

Singularities  of  the  flow  show  up  readily  with  the  method 
of  characteristics.   By  varying  the  shape  of  the  shock,  it  is 
possible  to  find  flows  with  either  an  envelope  of  characteris- 
tics in  the  supersonic  region  or  a  branch  point  in  the  subsonic 
region.   At  Mach  1.2,  singularities  always  occur  between  the 
shock  and  the  body,  while  at  Mach  1.5  a  flow  was  computed  with 
no  such  singularities. 
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1.   Introduction 

In  this  paper  the  flow  around  an  axially  symmetric  blunt 
body  in  a  steady  supersonic  flow  of  compressible  fluid  will  be 
discussed.  Our  aim  is  to  study  the  behavior  of  the  flow  between 
the  body  and  shock  for  low  Mach  numbers.  At  low  Mach  numbers 
the  detachment  distance,  i.e.,  the  distance  between  the  shock 
and  the  body,  is  relatively  large  and  hence  the  appearance  of 
singularities  in  the  flow  is  more  likely. 

It  is  physically  natural  to  prescribe  the  shape  of  the  body 
and  then  solve  the  equations  of  motion,  but  this  is  complicated 
by  the  fact  that  the  location  and  shape  of  the  shock  wave  deter- 
mined by  the  given  body  are  unknown.   Instead,  we  prescribe  the 
shock  and  find  the  generating  body  by  solving  a  Cauchy  problem 
for  the  equations  of  motion.   We  overcome  the  problem  that  the 
equations  of  motion  are  elliptic  near  the  vertex  of  the  shock 
wave  by  the  procedure  described  by  Garabedian  [3]>  which  exploits 
the  analyticity  of  the  problem  to  continue  the  equations  into 
the  complex  domain,  where  the  Cauchy  problem  becomes  properly 
posed. 

In  Section  2  we  present  the  mathematical  formulation  of 
the  problem  and  its  reduction  to  characteristic  normal  form. 
Section  3  deals  with  shock  conditions  and  the  choice  of  shock 
shape.   A  three-dimensional  method  of  attack  is  worked  out  in 
Section  4  with  a  discussion  of  the  corresponding  numerical 
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procedure.   We  present  in  Section  5  the  same  problem  solved  by 
the  method  of  characteristics  as  done  by  P.  Garabedian  [4]. 
The  superiority  of  this  method  at  low  Mach  numbers  becomes 
apparant.   By  plotting  the  characteristics,  the  location  of 
envelopes  or  "limiting  lines"  can  be  seen.   Also,  this  section 
contains  a  discussion  as  to  how  the  location  and  type  of  singu- 
larity depends  on  the  shock  shape  and  Mach  number. 

Much  work  has  been  done  on  the  two-dimensional  analog  of 
this  problem  by  Eva  V.  Swenson  [6,7].   Many  of  her  ideas  have 
been  utilized  in  this  paper.   In  the  axially  symmetric  case  we 
have  been  able  to  avoid  the  phenomenon  of  a  "hidden"  region  dis- 
cussed by  her,  and  to  achieve  flows  at  lower  Mach  numbers. 

The  author  is  deeply  grateful  for  the  assistance  of 
Professor  Paul  R.  Garabedian  who  suggested  the  problem,  and 
without  whom  this  paper  would  not  have  been  written. 
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2.   The  Equations  of  Motion  in  Normal  Form 

We  shall  start  with  the  equations  of  motion  for  steady 
axially  symmetric  flow  of  a  perfect,  compressible,  inviscid 
fluid  as  derived  in  [1].   Let  r  and  z  be  the  radial  and  axial 
coordinates  respectively,  and  let  u  and  v  be  the  velocities  in 
the  corresponding  directions.   Denoting  the  density  by  p  and 
the  pressure  by  p  we  have 

(1)  (rpu)z  +  (rpv)r  =  0 

which  expresses  the  fact  that  mass  is  conserved.   The  conser- 
vation of  momentum  is  expressed  by  the  two  equations 


(2a)  uuz  +  vur  +  i  pz  =  0 


(2b)  uvz  +  vvr  +  i  pr  =  0 


Equation  (l)  implies  the  existence  of  a  stream  function  1/ 
such  that 

(3)  ii^  =   rpu  ;        iz   =  -  rpv  . 

The   equation   of   state    is 

(*0  P  =  S(*)pr    ;  c2  =  §£=  ySWp7"1 
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where  y   is  a  gas  constant.   Conservation  of  energy  takes  the 
form  of  Bernoulli's  equation 

(5)  |(U2  +  v2)  +  ^E  =  h;   H  =  ^+^. 

From  equations    (l),    (2a, b),    (4)   and  (5),   we   are  able  to 
derive   the   system   (see   [1,6]) 

(6)  /c2-u2  0    \        /u\  /"-2uv  c2-v2>\      /u\     f£l.  -  ,VQ) 

0  c2-u2/        I  vi    +  l-(c2-u2)  0      J      \J+\-Q(c2-u2 


7(7-1)     S(f) 


Equations  (6)  and  (4)  comprise  the  equations  governing  the 
flow.  However,  for  our  purpose  it  will  be  convenient  to  intro- 
duce characteristic  coordinates.   Equation  (6)  is  in  the  form 

(7)  A  Uz  +  B  Ur  +  C  =  0 

and  hence  the  characteristic  directions  are  given  by  the  solution 
of 
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(8)  det  (AT  -  B)  =  0;      t  =  g 


Applying  this  to  (6),  we  find  that  there  are  two  roots  of  this 
equation,  namely 


,   A2   2   2 
-"UV  +  c/u  +v  -c 

(9)  \   "  ^^ 

—       c  -  u 


We  ignore  for  the  moment  the  possibility  that  the  discriminant, 

2   2   2 
u  +v  -c  ,  need  not  be  positive  and  will  return  to  this  point 

later. 

We  now  seek  to  find  a  vector  A  =  (A  ,A  )  fi   0  such  that 


(10)  A(At  -  B)  =  0 


This  leads  us  to  the  conclusion  that 


A„      2   2 

cm  £  =     c  ~v 

\±X'  A,    ,  2   2. 


vl    (c^-u^)t 

Multiplying  (7)  by  A  and  using  (10),  we  find  that 


(12)  *A[U_  +  tU  1  +  AC  =  0 


which  gives  us  two  scalar  equations  depend  1j  g  oi  which  valve  of 
t  is  used. 
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Let  ct(x,y)  =  const,  and  f3(x,y)  =  const,  be  the  characteris- 
tic ci  rves  with  slopes  x  and  t_  respectively.   Then  we  see 
from  the  definition  of  x  that 

(15)  ra   =  t+zq;      rp  =  x_zp 

so  that 


(14) 


U   =  z  (U  +  t,U  )  ; 
a    ax    z    +  r' 


UP  «  zp(Uz  +  T-Ur) 


Multiplying  (12)  by  z   and  letting  x  =  x  we  get,  with  the  aid 
of  (14), 


(15)  AAU   =  -z  AC  , 


and  similarly  for  x  =  x  we  see  that 


(16)  AAUg  =  -Zo^C 


We  now  insert  the  values  of  A  and  C  into  (15)  and  (l6)  and  use 
(ll)  with  A  =  x  to  derive 


2 

c  v 


(17a)   x+(c^)uq  +  (o*V>a  =  -ra[^  +  Q((c  -«  K   ~   uv)l 
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(17b)   T_(c2-u2)up  +  (c2-v2)vp  =  ~^[^-  +   Q((c2-u2)t_  -  uv)] 


Noting  that 


(18)       fa   =  ^zza  +  ^rra;      *p  =  ^zp  +  ^rp 


and  using  (4)  and  (13),  we  see  that 

(19)      |(^a  +  ^p)  =  rp[za(T+u-v)  -  zp(T_u-v)]  . 

Equations  (13),  (18)  and  (19)  comprise  the  characteristic 
normal  form  of  the  equations.   We  shall  return  to  them  later 
when  we  discuss  the  numerical  solution  of  this  system. 
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j5.   Initial  Conditions  and  Shock  Shape 

Since  in  our  problem  we  are  given  the  shock  shape,  it  is 
natural  for  us  to  consider  this  to  be  a  Cauchy  problem  with  the 
Cauchy  (or  initial  )  data  furnished  along  the  shock.   Data  along 
the  shock  can  be  calculated  once  the  free  stream  Mach  number  M„ 
and  the  shock  shape  z  =  f(r)  are  known.   The  derivations  are 
based  on  the  fact  that  mass,  momentum,  and  energy  are  conserved 
across  the  shock  (see  [1,6]). 

We  denote  values  in  front  of  the  shock  by  the  subscript  "0" 
and  values  behind  by  "1".   N  and  L  denote  the  normal  and  tangen- 
tial components  of  the  velocity  vector  with  respect  to  the  shock, 
Choosing  a  scale  so  that  p„=  p0  =  1,  we  find  that 

Ll  =  L0  ' 
(21)  Pl  -  N2(l  -  ix2)    -  u2;    |x2=  *=1  , 

Pi  =  VN1  • 

We  shall  take  the  direction  of  incident  flow  to  be  parallel  to 
the  axis  of  the  body,  which  is  the  z-axis.   Hence  vQ  =  0  and 
u„  =  Un-   Since  we  have  chosen  p„  =  p„  =  1  we  get  from  (5)  that 


(22)  c0  =  T 


The  Mach  number  M  is  defined  as  the  ratio  of  the  spe-ed  of 
the  fluid  to  the  speed  of  sound.   That  is, 


/2,  2 

(23)  M=^L_ 


and  hence  from  (22)  we  see  that  UQ  =/y'M0. 

Let  9   be  the  angle  which  the  normal  to  the  shock  wave  makes 
with  the  horizontal.   Then  for  the  shock  wave  z  =  f(r)  we  have 

(24)  LQ  -  UQsin  9    ;      NQ  =  UQcos  0 


where 


(24")    sin  9   ~   f'(r)     ;  cos  9   -     X               . 

/l  +  f ,2(r)  /l+f ,2(r) 

Using  (2l)-(24')  we  can  obtain  the  following  expressions  for 
1^,  Nlf  p1  and  p1: 


(25) 


Lx  =  UQ   sin  9    , 


p1   =    (1-|x2)UqCos20   -  \?    , 


Nl   =   [f+T  +   M-2U2cos28]  /  U0cos   9    , 
Pl  =  u2cos20[^+   n2U2cos2e]-1    . 
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Finally  we  obtain  vu  and  v-.  from  the  relations 


(26) 


u.,  =  N-.  cos  6  +   L-,  sin  0 


v-,  =  L-,  cos  9   -   N-,  sin  9    . 


2     2    2 
As  remarked  earlier,  the  discriminant  u  +  v   -  c   need 

not  be  positive,  in  which  case  the  characteristics  are  not  real. 
We  see  from  (23)  that  the  discriminant  is  positive,  zero,  or 
negative  according  to  whether  M>  1,  M=l,  or  M  <  1.   In  the 
case  that  M  <  1,  the  Cauchy  problem  is  not  well  posed.   However, 
by  our  assumption  of  analyticity  of  the  shock  shape  and  the 
elementary  from  of  (25)  and  (26),  we  can  consider  r  as  a  complex 
variable  and  extend  the  Cauchy  data  into  the  complex  domain. 
Through  every  point  in  complex  (r,z)  space  there  are  two  dis- 
tinct complex  characteristics  if  M  £   1.   Since,  as  we  shall  see, 
our  method  of  solution  only  requires  that  the  two  characteristics 
be  distinct,  we  find  this  method  of  complex  extension  [3,^,6] 
extremely  useful.   The  points  where  M  =  1  will  be  of  special 
interest.   The  set  of  real  points  for  which  M  =  1  will  be  called 
the  sonic  line. 

We  are  particularly  interested  in  computing  flows  at  low 
Mach  numbers.   To  have  any  physical  significance  the  flow  must 
not  have  any  singularities  in  the  region  between  the  shock  and 
the  body,  which  is  defined  by  the  streamline  ii   =  0.   When  this 
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is  achieved  we  say  that  we  have  "smooth"  flow. 

In  the  selection  of  the  shock  shape  z  =  f(r)  we  must  satisfy 
the  following  requirements: 

(a)  f (r)  has  to  be  an  even  function; 

(b)  f(r)  has  to  be  everywhere  non-characteristic; 

(c)  f (r)  has  to  be  analytic  in  a  relatively  large  region 
about  the  real  r  axis,  and  real  on  the  real  axis; 

(d)  f(0)  =  0. 

Since  the  characteristics  have  shapes  t+  =  +  /Mq-1  asymptotically, 
and  since  (a)  implies  that  f'(0)  =  0,  we  can  conclude  from  (b) 
that  |f'(r)|  <  /m^-1  for  all  real  r. 

Since  we  want  to  find  the  solution  in  a  large  region,  we 
do  not  want  any  singularities  close  to  the  real  axis.  The 
function 

(27)  z  =  f (r)  =  a(/p2+r2  -  p) 

is  perhaps  one  of  the  simplest  of  such  functions.   Taking  the 
derivative  we  discover  that 

(28)  f'(r)  =   ar(p2  +  r2)"l/2 


from  which  we  see  that  lim  f ' (r)  =  a.   Also,  by  computing  the 

r->oo 
second  derivative,  we  note  that 
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(29)  f"(r)  =  ap2(p2  +  r2)~5/2  , 

and  hence  the  curvature  at  r  =  0  is  a/p • 

In  our  work  we  found  it  convenient  to  let 

(30)  a  =  A  /m2-1  ;      (3  =  AMQ 
in  which  case  (27)  becomes 

(31)  z  =  f(r)  =  A/M2-l(/(AM0)2+r2   -  AMQ) 

-1   2    l/2 
which  has  curvature  MQ  (MQ-l)  '       at  the  origin  independent  of 

the  value  of  A.   For  0  <  A  <  1  it  also  has  the  property  that 

the  curve  is  everywhere  non-characteristic. 
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4.   A  Three -Dimensional  Method 

Instead  of  transforming  into  characteristic  coordinates, 
it  is  possible  to  attack  the  problem  directly.   The  system  (6) 
along  with  (3)  can  be  written  as 

2      2  X  x  2 

2uv  _  c    -v  0  \      .  .  XuvQ  _  c  v 

'  c2-u2  \f\        /c2-u2  r(c2-u2) 


(32)|    v   |-    l  0  0  Q 

-  rpv 

which  has  the  form 

(33)  Uz  =  RUr  +  S  . 

Letting  £  =  r  and  £  =  z-f (r)  we  can  transform  the  shock 
into  a  coordinate  axis.   Equation  (33)  now  takes  the  form 


(3^)      U.  -  [I  +  f«(e)R]_1RU^  +  [I  +  f'C^R]"^ 


and  a  simple   calculation  yields 

rc2-u2  f(£)(c2-v2)  0 

(35)      [I  +  f (l)R]"1  -  J  \    -f'(^)(c2-u2)  (c2-u2)  +  2uvf'(!)  0 

0  0  D 
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2   2. 


2   2N„t2 


with  D  =  (c  -u  )  +  2uvf'(£)  +  (c  -v  )f '*(£) .   System  (52)  now 
becomes 


(36) 


/u\ 


/2uv+(c2-v2)f  •    -(c2-v2) 


0 


V 


D 


2   2 
c  -u 


(c2~v2)f     0 


U 


0 


0  0 

Q[uv+(c2-v2)f]  -  41- 
+  J  |  Q[uvf <+c2-u2]  +  ifL 
-   IpvD 

It  can  be  shown  that  (36)  is  not  a  hyperbolic  partial 

2    2    2 
differential  equation  when  u  +  v   -  c   is  negative.   In  order 

to  make  the  system  symmetric  hyperbolic,  we  perform  a  complex 

extension  of  the  ^.-coordinate ,  i.e.,  the  |-axis  becomes  the 

real  axis  of  a  complex  plane  perpendicular  to  the  £  axis. 

Under  this  complex  extension  the  system  (36),  which  is  of  the 

form 


(37) 


U.  =  RUe  +  C  , 


becomes 


-1^- 


R+R  ntt   .  ,R-R 


(38)  U^  =  (^-)ue  +  (^-)UT]  + 


or 


(39)  U^  =  AU^  +  BU^  +  C 


where  £,  is  now  replaced  by  £  +  ±t\   in  (36). 

The  above  procedure  has  the  advantage  of  making  the  problem 
well  posed,  since  A  and  B  are  both  Hermitian;  however,  it  also 
increases  the  number  of  independent  variables  from  two  to  three. 
Thus,  we  have  changed  our  two-dimensional  problem  into  a  three- 
dimensional  one. 

For  the  numerical  solution  of  this  three-dimensional  prob- 
lem, the  Lax-Wendroff  method  is  used  [2].   In  this  case,  £  plays 
the  role  of  the  time  variable. 

The  stability  criterion  for  (39)  is  arrived  at  by  assuming 
that  A  and  B  have  constant  coefficients.   By  substituting 
V  *=  VQ  einA*+:LinAn  into  the  Lax-Wendroff  difference  equations, 
we  get  the  "symbol,"  or  amplification  matrix, 


(40)   Q(A1,A2,a,p)  =  I  +  i(cos  a   cos  P)(A-LA  sin  a  +  AgB  sin  6) 


-  2[AXA  sin  a  +  A2B  sin  p]2  +  A£C 


where 
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a  =  nA£  ;        (3  =  mAr\ 


and 


>    AC       -x    Ar 

1    AT  '         2    Ati 


Let 


(41)  P  =  A  A  cos  0  +  A2B  cos  9 


where 


cos  9   =    sln  a    ■  ;      sin  0  =    sln  P 

i 2 2~~  '  / 2" 2~ 

/sin  a+sin  p  /sin  a+sin  f 


We  can  now  rewrite  (40)  as 


a  / 5 p~~ 

(42)  Q  =   I  +   i(cos  a  +   cos  p)  /sin  a+sin  p   F 

-   2(sin2a  +    sin2p)F2   +   A£C    . 


A 

We  denote  the  eigenvalues  of  Q  by  q.  and  those  of  F  by 
[i. ,    i  =  1,2,3-   A  necessary  condition  for  stability  is  that 
|q.|  <  1  +  0(A£)>  i  =  1,2,3-   Noting  that  the  term  A£  C  is 
0(A£)  and  diagonalizing  (42),  we  arrive  at  the  inequality 
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,  2Qx2  ,.2r«..2 


(43)      |q1T  <  1  +  2(sin^a  +  sin^p)^  M^[2^  -  1] 


.2 


from  which  we  see  that  |q.  |  <  1  only  if  p..  <  l/2.   By  recalling 
(4l)  we  see  that  a  necessary  condition  for  stability  is  that 
the  spectral  radius 


(44)         p[(\&   cos  6  +   *2B  sin  0)2]  <  i 


for  all  0.  By  a  simple  calculation  we  can  show  that  this  will 
indeed  be  the  case  if 


(*5)  ^||a||2  +  A2||B||2  <\ 


where  the  norm  used  here  is  the  spectral  norm. 

A  program  has  been  written  for  the  CDC  6600  which  solves 
the  problem  by  this  method.   Equation  (31)  was  used  as  the 
shock  shape  with  A  =  1.   Initially  a  triangular  grid  of  data, 
with  r\   >  0,  is  prescribed  as  shown  in  Figure  1.   Points  for 
negative  values  of  r\   are  obtained  by  the  Schwarz  reflection 
principle.   The  stability  criterion  is  checked  at  every  point 
on  the  grid,  and  the  £-step  is  chosen  small  enough  to  satisfy 
it  everywhere. 

This  method  worked  well  for  Mach  numbers  MQ  near  6,  but 
when  the  Mach  number  was  significantly  lowered,  it  was  not 
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possible  to  satisfy  the  stability  condition  and  still  reach 
the  body.   Since  the  stability  criterion  was  more  easily  satis- 
fied at  some  points  in  the  grid  than  others,  a  substitution  was 
made  in  the  independent  variables  in  hopes  of  reaching  the  body 
with  a  stable  run.   The  transformations  made  were  of  the  form 


(46) 


£'  =  a^3  +  b|  , 


3 

r)«  =  ct}v  +  drj 


with  various  values  of  a,b,c,  and  d.  Under  the  transformation 
equation  (39)  takes  the  form 


U*  ,    +  B  - 

d£'/d£    ^  dTi'/dn       ' 


(47)  U^  =  A    X      ^  ,  +  B i—  Un ,  +  C 


so  that  in  effect  A  and  B  are  replaced  by 


-| and  1- 


3a£,  +  b       3ct}  +  d 

Much  of  our  work  was  done  at  Mach  2.0.   Many  runs  were 
made  with  various  values  of  a,b,c  and  d  until  we  were  able  to 
compute  the  whole  subsonic  portion  of  the  flow  without  violating 
our  stability  condition.   The  reason  for  this  difficulty  became 
apparant ,  however,  when  the  method  of  characteristics  (to  be 
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explained  later)  exhibited  a  limiting  line  just  past  the  region 
in  which  we  were  able  to  solve. 

A  typical  run  on  the  CDC  6600  with  97  points  initially  on 
the  real  £,  axis  takes  about  10  minutes  of  CP  time.   The  program 
uses  only  core  storage  and  requires  about  50,000  words  of  memory. 
The  machine  carries  about  14  digits  for  all  calculations  so  it 
is  assumed  that  round-off  error  is  insignificant.   Results  of  a 
flow  at  Mach  2.0  were  tabulated.   Table  1  shows  the  convergence 
at  a  particular  point,  the  last  line  being  the  significant 
figures  after  extrapolation.   Finally  this  flow  is  plotted  in 
Figure  2  showing  the  shock,  the  body,  the  sonic  line,  the  flow 
direction  and  the  limiting  line. 
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5-   The  Method  of  Characteristics 

In  Section  2  the  equations  of  motion  were  put  into  charac- 
teristic normal  form.   We  will  now  discuss  the  numerical  solu- 
tion of  these  equations  by  the  method  of  Massau.   The  basic 
idea  of  the  method  is  that  given  any  two  points  in  r,z  space, 
and  the  values  u,v  and  if   there,  we  can  compute  a  new  point  and 
the  solution  at  that  point  to  second  order  accuracy  by  looking 
at  the  a  and  p  characteristics  passing  through  the  two  points. 
A  more  detailed  analysis  of  the  method  of  Massau  can  be  found 
in  [2,6].   We  should  also  point  out  that  this  method  works  pro- 
viding the  two  characteristics  intersect.  This  will  indeed  be 
the  case  provided  the  a  and  (3  characteristics  are  not  parallel. 

Since  this  method  is  applied  on  neighboring  points  in  a  grid, 

2    2    2 
we  see  that  this  will  happen  only  when  u  +v  =  c  ,  or  M  =  1. 

By  allowing  for  the  possibility  that  the  characteristics  may 
be  complex  we  circumvent  the  difficulty  of  the  equation  being 
elliptic  near  the  origin.  Thus  our  only  difficulty  in  computa- 
tion will  occur  at  points  where  M  =  1. 

In  our  problem  we  have  extended  the  initial  data  into  the 
complex  plane  by  analytic  continuation.   The  form  of  the  shock 
shape  (31)  and  the  elementary  form  of  (25)  and  (26)  gives  us 
an  explicit  expression  for  the  initial  data  throughout  the 
complex  plane.   Along  that  portion  of  the  real  shock  where 
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M  >  1  the  characteristics  will  be  real  and  distinct,  while  the 
portion  where  M  <  1  has  characteristics  occurring  as  complex 
conjugates.   Data  is  prescribed  along  a  chosen  path  in  the  com- 
plex r-plane .   There  is  no  restriction  on  the  path  except  that 
it  be  continuous.   This  determines  a  wedge  of  data  in  complex 
r-z  space  bounded  by  an  a  characteristic  passing  through  one 
end  of  the  initial  path  and  a  B  characteristic  passing  through 
the  other  end  (see  Figure  3) . 

For  the  physical  problem  we  are  only  concerned  with  finding 
solutions  at  points  where  r  and  z  are  real.   The  fact  that  the 
characteristics  are  real  along  the  real  r-axis  where  M  >  1 
implies  that,  by  choosing  an  initial  path  on  this  line,  the 
whole  domain  of  solution  will  be  real  (see  Figure  4).   We  also 
note  that,  for  the  subsonic  portion,  the  a  and  P  characteristics 
through  any  conjugate  points  are  complex  conjugates  and  hence 
must  intersect  at  a  real  point.   By  choosing  a  symmetric  path 
about  the  real  r-axis  in  the  subsonic  region  we  obtain  one  line 
of  real  data  within  our  wedge  of  complex  solution  (see  Figure  5)« 
Referring  to  Figure  5,  we  see  that  the  two  characteristics 
passing  through  point  P  emanate  from  a  and  a.      Similarly,  through 
any  point  with  M  £   1  there  are  two  characteristics.   To  reach 
any  point  Q,  all  we  need  do  is  trace  the  characteristics  through 
Q,  back  to  the  initial  plane  and  apply  our  method  to  some  initial 
path  connecting  them.   The  initial  path  must,  of  course,  have 
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the  property  that  the  wedge  of  data  does  not  contain  points  for 
which  M  =  1. 

There  is  another  possibility,  however.   In  tracing  the 
characteristics  back  to  the  initial  r-plane  one  of  the  charac- 
teristics may  pass  through  a  point  where  M  -  1.   In  this  case 
there  is  said  to  be  a  "hidden"  region  and  all  points  with  this 
property  are  said  to  be  in  the  hidden  region.   For  the  two 
dimensional  analog  of  this  problem,  this  phenomenon  cannot  be 
avoided  by  any  choice  of  the  shock  shape.   A  procedure  to  over- 
come this  difficulty  and  to  get  into  the  hidden  region  has  been 
worked  out  by  E.  Swenson  [7].   Application  of  this  procedure 
has  been  successful,  but  by  clever  choice  of  the  shock  shape  it 
is  possible  to  avoid  the  phenomenon  of  a  hidden  region  in  the 
axially  symmetric  case. 

We  now  turn  to  another  aspect  of  the  procedure.   We  must 
also  consider  the  prospect  of  two  characteristics  of  the  same 
family  intersecting.   If  they  do  then  the  solution  is  not  single 
valued  at  the  point  of  Intersection.   In  the  supersonic  region 
the  characteristics  might  envelope,  forming  a  limiting  line  for 
the  solution.   In  fact,  at  low  Mach  numbers  this  difficulty 
becomes  paramount  to  solving  the  problem,  since  the  increased 
detachment  distance  makes  the  likelihood  of  singularities 
greater.  The  remainder  of  this  paper  will  be  devoted  to  a 
discussion  of  such  singular  behavior  and  its  dependence  on  the 
shape  of  the  shock. 
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We  did  a  great  deal  of  experimentation  at  various  Mach 
numbers  using  (31)  as  the  shock  shape  with  various  values  of  A. 
We  first  consider  the  case  where  A  =  1.   For  the  three-dimensional 
method  this  was  the  only  value  of  A  which  was  tried.   From  (28) 
and  (30)  we  see  that  the  shock  becomes  characteristic  asympto- 
tically.  The  solution  always  contains  an  envelope  of  charac- 
teristics extending  to  infinity.   This  envelope  is  a  limiting 
line  for  the  solution.   That  is,  a  single-valued  solution  exists 
only  as  far  as  the  limiting  line.   Depending  on  the  Mach  number, 
the  other  end  of  the  envelope  either  has  to  cross  the  body  be- 
fore it  meets  the  sonic  line  or  meet  the  sonic  line  first. 
Because  of  the  increased  detachment  distance  at  low  Mach  numbers, 
the  latter  is  the  case  and  hence  the  complete  body  cannot  be 
found.   At  Mach  1.5  this  turned  out  to  be  the  case.   At  Mach  2.0, 
we  see  in  Figure  2  that  the  limiting  line  reaches  the  body 
first  and  occurs  just  beyond  the  region  we  were  able  to  solve 
for  with  the  three-dimensional  method.   This  explains  the  dif- 
ficulty in  trying  to  increase  the  size  of  the  region  of  solution 
by  a  change  of  variables  using  the  three-dimensional  method. 

It  was  hoped  that  this  difficulty  would  be  eliminated  by 
lowering  the  value  of  A.   When  the  value  of  A  was  lowered 
slightly,  the  envelope  no  longer  extended  to  infinity.   The 
length  of  the  singularity  decreased  as  A  was  lowered  still  more. 
The  place  where  the  other  end  meets  the  sonic  line  does  not 
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change  significantly.   Figure  6a-6c  shows  three  flows  with 
different  values  of  A  for  Mach  1.5.   As  A  is  decreased  the 
length  of  the  limiting  line  decreases,  but  the  other  end  still 
meets  the  sonic  line  in  front  of  the  body. 

As  the  value  of  A  is  decreased  more,  this  envelope  com- 
pletely disappears.   If  there  is  a  hidden  region  in  the  flow, 
the  envelope  may  be  in  it,  however,  so  that  we  prefer  to  find 
flows  without  hidden  regions.   With  Mach  1.5  and  A  =  .75  we 
were  able  to  compute  a  flow  with  a  hidden  region.   By  a  method 
due  to  Swenson  [7],   we  were  able  to  show  that  there  was  no 
limiting  line  inside  the  hidden  region.   At  A  =  -70  the  hidden 
region  completely  disappeared  and  we  were  able  to  compute  a 
smooth  flow.   In  this  case  the  sonic  line  is  in  two  pieces, 
each  of  which  extends  to  infinity.  The  calculations  were  done 
on  the  CDC  6600  with  two  programs,  one  for  the  supersonic  region 
and  one  for  the  subsonic  region.   The  former  uses  real  arithmetic 
and  takes  about  80  seconds  CP  time  for  an  initial  grid  of  321 
points,  while  the  latter  uses  complex  arithmetic  and  takes  500 
seconds  CP  time  for  the  same  sized  grid.   Several  runs  of  the 
second  program  were  needed,  since  each  run  yields  only  one  line 
of  solution.   Away  from  the  sonic  line,  five  or  six  significant 
figures  were  achieved  while  points  on  the  sonic  line  were  extra- 
polated to  only  three  figures.   Table  2a  shows  convergence  at  a 
point  in  the  supersonic  region,  while  Table  2b  shows  convergence 
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at  a  point  in  the  subsonic  region.   The  last  line  in  each  table 
is  obtained  by  extrapolation.   Finally,  Figure  7  is  a  plot  of 
this  flow  displaying  the  shock,  body,  sonic  lines,  and  the 
characteristics  in  the  supersonic  region. 

At  Mach  1.2  the  value  of  A  was  lowered  until  the  envelope 
disappeared.   In  this  case  these  were  singularities  in  the 
hidden  region.   When  the  value  of  A  was  lowered  more,  a  branch 
point  developed  in  the  subsonic  region  in  front  of  the  body 
(see  Figures  8a,  8b).   As  A  is  decreased  to  zero,  the  branch 
point  moves  closer  to  the  z-axis,  always  remaining  in  front  of 
the  body. 
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Table   1.       (Convergence) 
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Table  2A.   (Convergence) 
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Table  2B.   (Convergence) 
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Fig.    1.      Typical   grid  used   with  7  points   on   the   initial   real   line 
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Fig.  2.   Flow  at  Mach  2.0  computed  by  the  three-dimensional  method 
with  7=1.4,  A  =  1.00 
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Fig.  5.   The  domain  of  solution  of  a  given  initial  path 
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Fig.  4.   An  initial  path  on  the  real  axis  in  the  supersonic  region  and 
the  corresponding  domain  of  solution  determined  by  this  data 
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Fig.  5.   A  symmetric  path  about  a  point  r  on  the  real  r  axis  in  the 
subsonic  region  and  the  line  of  solution  in  the  real  r,z 
plane  which  it  yields 
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Fig.  6a 


Figs.  6a-6c, 


Three  flows  at  Mach  1.5,  7  =  1.4,  with  different  values 
of  A.   The  envelope   of  characteristics  forms  a  limiting 


line, 
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Fie.    6b 
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Fig.    6c 
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Fig.  7.   Flow  computed  at  Mach  1.5  with  7  =  1.4  and  A  =  .70 
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Fig.  8a 


Figs.  8a-8b, 


Initial  path  and  path  of  real  solution  which  indicates  a 
branch  point.   Flow  Is  at  Mach  1.2,  y   =  1.4,  and  A  =  . 40 
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